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FOREWORD 


Opinions,  interpretations,  conclusions  and  recommendations  are 
those  of  the  author  and  are  not  necessarily  endorsed  by  the  U.S 

Army. 

Where  copyrighted  material  is  quoted,  permission  has  been 
obtained  to  use  such  material. 


Where  material  from  documents  designated  for  limited 
distribution  is  quoted,  permission  has  been  obtained  to  use  the 

material . 

Citations  of  commercial  organizations  and  trade  names  in 
this  report  do  not  constitute  an  official  Department  of  Army 
endorsement  or  approval  of  the  products  or  services  of  these 
organizations . 

In  conducting  research  using  animals,  the  investigator  (s) 
idhered  to  the  "Guide  for  the  Care  and  Use  of  Laboratory 
Animals,"  prepared  by  the  Committee  on  Care  and  use  of  Laboratory 
Animals  of  the  Institute  of  Laboratory  Resources,  National 
Research  Council  (NIH  Publication  No.  86-23,  Revised  1985) . 

For  the  protection  of  human  subjects,  the  investigator (s) 
.idhered  to  policies  of  applicable  Federal  Law  45  CFR  46. 

In  conducting  research  utilizing  recombinant  DNA  technolop, 
the”investigator (s)  adhered  to  current  guidelines  promulgated  by 
the  National  Institutes  of  Health. 

In  the  conduct  of  research  utilizing  recombinant  DNA,  the 
investigator (s)  adhered  to  the  NIH  Guidelines  for  Research 
Involving  Recombinant  DNA  Molecules . 

In  the  conduct  of  research  involving  hazardous  organisms, 
the"investigator (s)  adhered  to  the  CDC-NIH  Guide  for  Biosafety  in 
Microbiological  and  Biomedical  Laboratories. 
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Introduction 

The  military  employs  a  broad  range  of  laser  radiation  in  training  devices, 
rangefinders,  target  designators,  communications  devices  and  other  instruments. 
This  equipment  emits  either  single  pulses  or  sequences  of  pulses  in  beams  of 
various  diameters.  The  research  performed  under  this  contract  directly  supports 
the  U.  S.  Army  Medical  Research  and  Materiel  Command  (USAMRMC)  mission 
to  assess  the  health  effects  and  hazards  of  nonionizing  electromagnetic  radiation 
from  such  laser  systems.  The  data  obtained  will  support  evaluation  of  current 
permissible  exposure  limits  and  will  aid  health  pohcy  makers,  both  within  and 
outside  the  DoD,  in  developing  injury  prevention  criteria.  The  general  approach  is 
to  make  direct  determinations  of  damage  threshold  levels  for  non-ionizing 
radiation  for  specific  exposure  conditions  (e.g.,  wavelengths,  pulse  durations,  etc.) 
and  to  develop  models  of  the  damage  mechanism  that  enable  the  extension  of  the 
results  to  other  exposure  conditions. 

Under  past  support  from  the  Army  Medical  Research  and  Development 
Command  we  determined  corneal  damage  thresholds  for  CO2  lasers  that  emit 
single  and  sequences  of  pulses  having  durations  of  1  ms  and  longer  and  developed 
thermal  damage  models  for  predicting  this  type  of  damagei-^.  We  also 
determined  corneal  damage  thresholds  for  single  80  ns  pulses  of  CO2  laser 
radiation  and  for  three  multiple-pulse  (2  pulses  at  1  Hz.;  and  2  and  8  pulses  at  10 
Hz.).6>  7  In  the  case  of  these  very  short  pulses,  damage  mechanisms  other  than 
thermal  (e.g.,  acoustical  pressure  pulses)  could  also  play  a  role.  Light  and 
electron  microscopy  revealed  unusual  disruptions  of  the  anterior  epithelial 
surface  for  the  threshold  single-pulse  exposure.  The  characteristics  of  these 
disruptions  differed  from  those  observed  with  simple  thermal  damage  at  longer 
pulse  durations  and  appeared  to  be  consistent  with  a  mechanical  [e.g.,  acoustic] 
damage  mechanism.  However,  the  calculated  temperature  increase  produced  by 
the  threshold  exposure  was  only  slightly  lower  than  that  calculated  for  threshold 
exposures  having  durations  of  1  ms  and  longer.^*  ^  Thus  we  could  not  exclude  a 
thermal  damage  mechanism,  with  the  sharp  temperature  gradients  leading  to 
marked  differences  in  the  character  of  the  damage  as  compared  to  damage  from 
longer  duration  exposures.  In  the  case  of  the  preliminary  multiple-pulse  data 
noted  above,  the  10  Hz  results  contained  an  (apparent)  discrepancy  in  that  the 
threshold  for  the  8  pulse  sequence  was  higher  than  that  for  the  2  pulse  sequence. 

This  study  aims  to  expand  the  database  of  thresholds  for  sequences  of  80  ns 
of  CO2  laser  pulses  and  to  refine  the  preliminary  10  Hz.  thresholds.  In  addition  it 
will  attempt  to  clarify  the  damage  mechanism(s)  for  these  t5q)es  of  pulses  and 
ultimately  will  begin  to  address  damage  from  mid-infrared  wavelen^hs  where 
the  radiation  is  more  penetrating. 


Methods 

Short-pulse  exposures  are  made  with  a  Boston  Laser  (Model  220S)  CO2-TEA 
laser  operated  in  the  TEMoo  mode.  This  laser  delivers  80  ns  pulses  at  pulse 
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repetition  frequencies  up  to  20  Hz.  Mode  quality  is  verified  and  the  beam  diameter 
is  measured  at  the  beginning  of  each  experimental  session  using  a  Spiricon 
linear  pyroelectric  array.  The  detector  has  64  elements  on  200  |jm  centers.  It  is 
mounted  on  a  vertical  micropositioner  and  is  read  out  with  a  LeCroy  9354M  digital 

oscilloscope.  These  latter  two  pieces 
of  equipment  are  major 
improvements  to  the  apparatus 
since  our  previous  contract.  Figure 
1  shows  the  profile  of  a  beam  having 
a  1/e  diameter  of  3.50  mm.  Pulse 
energy  is  measured  with  a 
Scientech  detector  immediately 
before  and  immediately  after  each 
exposure. 


New  Zealand  white  rabbits  of 
either  sex  weighing  4-5  pounds  are 
used  for  the  experiments.  The 
rabbits  are  treated  in  accordance 
with  the  Guide  for  the  Care  and  Use 
of  Laboratory  Animals  (DHEW 

Figure  1 :  Profile  of  TEMqo  beam  d^e  =  3.50  mm.  Publication  No.  (NIH)  85-23,  Revised 

Edition,  1985  and  with  the 
Association  for  Research  in  Vision 
and  Ophthalmology  Resolution  on  the  Use  of  Animals  in  Research.  Prior  to 
exposure  the  rabbits  are  anesthetized  with  an  intramuscular  injection  of  xylazine 
and  ketamine  hydrochloride  (Rompun-Ketaset)  in  the  proportions:  60%  of  20 
mg/ml  Rompun  to  40%  of  100  mg/ml  Ketaset  by  volume.  In  addition  a  topical 
anesthesia  (tetracaine  hydrochloride  1/2%)  is  applied  to  each  eye  before  exposure 
-  and  a  drop  of  homatropine  bromide  5%  is  instilled  to  dilate  the  pupil.  This 
facilitates  examination  of  the  exposed  corneas  for  minimal  lesions.  The 
anesthetized  animals  are  placed  in  a  conventional  holder  for  exposure.  A 
speculum  is  inserted  in  the  eye  about  one  minute  before  exposure  and  the  eye 
irrigated  with  BSS  solution  (Alcon  Surgical);  however,  in  order  to  create  a 
reproducible  “tear  film,”  the  irrigation  is  stopped  20  sec  before  the  exposure  and 
excess  fluid  is  blotted  at  the  limbus.  One  exposure  is  made  to  each  eye.  One-half 
hour  after  exposure  the  rabbits,  still  under  anesthesia,  are  sacrificed  with 
Beuthanasia-D  administered  in  an  ear  vein.  The  eyes  are  enucleated  and 
examined  for  damage  using  a  Nikon  FS-3  photo  slit-lamp. 

The  criterion  that  we  use  for  minimal  epithelial  damage  is  that  due  to 
Brownell  and  Stuck^,  namely  the  presence  of  a  superficial  gray-white  spot  that 
develops  within  one-half  hour  after  exposure.  We  have  found  that  the  damage 
threshold  is  sharply  defined;  i.e.,  only  rarely  is  there  overlap  between  exposures 
that  produce  minimal  lesions  and  those  that  do  not.  Therefore  we  do  not  use 
statistical  procedures  such  as  probit  analysis  in  order  to  determine  the  threshold, 
as  these  would  require  the  use  of  more  animals  than  we  deem  necessary.  We 
make  one  exposure  per  eye,  bracketing  exposures  above  and  below  threshold.  The 
bracket  is  narrowed  until  there  is  only  about  a  10%  difference  in  energy  between 
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an  exposure  that  produces  a  minimal  lesion  and  one  that  does  not.  The  threshold 
is  taken  to  be  at  the  center  of  the  bracket. 

Temperature  calculations  are  based  on  a  Green’s  function  solution  to  the 
heat  conduction  equation  for  an  incident  beam  with  a  Gaussian  irradiance  profile 
that  is  absorbed  according  to  Beer’s  law.  The  beam  is  assumed  to  impinge  on  a 
semi-infinite  slab  and' to’ have  a  constant,  peak  irradiance  for  the  duration  of  the 
exposure.  We  also  assume  that  conduction  is  the  only  mode  of  heat  transfer  and 
that  no  heat  is  lost  to  the  air  at  the  epithelial  interface.  The  solution  T{r,z,t),  where 
r  is  the  radial  distance  from  the  beam  axis,  z  is  the  depth  into  the  cornea,  and  t  is 
time, has  the  form  of  a  definite  integral  that  can  be  evaluated  numerically.^’  ^ 
We  wrote  a  program  to  evaluate  this  integral  under  our  initial  contract.  The 
program  ran  on  the  APL  mainframe  computer  which  is  no  longer  in  operation. 
We  have  modified  the  original  program  as  discussed  below  so  that  it  runs  on  a 
desktop  computer  (in  our  case  a  Power  Macintosh).  The  modifications  are 
discussed  below  and  a  copy  of  the  new  program  is  listed  in  Appendix  2. 

Results  and  Discussion 

The  contract  with  USAMRMC  was  signed  November  9,  1995.  Budgets  were 
established  at  APL  by  the  middle  of  December.  During  January  we  (Drs.  R.  L. 
McCally  and  C.  B.  Bargeron)  conferred  with  the  Contracting  Officer’s 
Representative,  Mr.  Bruce  Stuck,  to  reaffirm  his  research  needs  and  priorities.  In 
this  regard  it  should  be  noted  that  the  Statement  of  Work  for  the  first  year  dated 
February  24, 1995  was  the  same  as  in  the  original  1991  proposal  and  we  wanted  to 
confirm  that  the  research  described  there  was  still  appropriate.  That  Statement  of 
Work  was  as  follows: 

1.  The  data  base  of  threshold  conditions  for  80  ns  pulses  from  CO2-TEA 
lasers  will  be  extended  by  measuring  the  threshold  energy  densities  for 
sequences  of  32,  128  and  1024  pulses  at  10  Hz  and  sequences  of  2,  8,  32, 
and  128  pulses  at  20  Hz.  The  existing  thresholds  for  2  and  8  pulses  at 
10  Hz  will  be  refined. 

2.  An  understanding  of  the  damage  mechanism  for  such  pulses  will  be 
pursued  by: 

a)  determining  how  lowering  the  temperature  of  the  epithelium 
affects  the  threshold, 

b)  obtaining  light  and  electron  micrographs  of  the  damage, 

c)  using  high-speed  photography  to  investigate  ablation  of  material 
from  the  corneal  surface,  and 

d)  measuring  pressure  transients  in  model  systems. 

3.  The  required  characteristics  of  potential  near-infrared  laser  source(s) 
for  damage  studies  will  be  identified  by  performing  detailed 
temperature  calculations  in  the  1.3  to  2.5  (im  wavelength  region. 
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Mr.  Stuck  confirmed  that  the  multiple-pulse  damage  threshold 
determinations  noted  there  were  still  of  major  interest.  Because  our  laboratory 
had  not  been  used  for  this  type  of  research  since  termination  of  our  last  contract 
in  1989  many  of  the  pieces  of  equipment  had  migrated  to  other  experiments.  Thus 
we  began  a  major  effort  to  locate  them  and  re-assemble  the  apparatus  so  that  the 
threshold  determinations  could  be  made. 

As  noted  above  in  the  Methods  Section,  the  program  for  calculating 
temperature  histories  that  we  had  developed  during  our  earlier  contract  (1977-89) 
had  to  be  revised  because  the  mainframe  computer  on  which  it  was  designed  to 
run  is  no  longer  in  operation.  The  original  program  used  mainframe  library 
routines  for  calculating  complementary  error  functions  and  the  routine  DCADRE 
for  evaluating  the  integrals  that  give  temperature  as  a  function  of  time  and 
position  in  the  cornea.  The  routine  DCADRE  was  an  adaptive  Gaussian 
quadrature  method  for  approximating  integrals  that  also  was  no  longer  available. 
In  revising  the  program  for  the  Power  Macintosh  using  FORTRAN  77  (LS 
FORTRAN,  Fortner  Research  LLC)  we  used  routines  obtained  from  Numerical 
Recipes  *0  for  evaluating  the  complementary  error  functions  and  the  IMSL  routine 
DQDAG  for  evaluating  the  integrals.  The  routine  DQDAG  uses  a  globally 
adaptive  scheme  that  is  based  on  Gauss-Konrad  rules.  We  checked  the  accuracy 
of  the  error  function  calculations  by  comparing  the  results  to  values  computed 
using  Mathematica  (Wolfram  Research).  The  results  agreed  to  double  precision 
accuracy  (15  significant  figures).  The  entire  program  was  tested  by  repeating 
several  calculations  done  previously  on  the  mainframe.  Tests  were  done  for  both 
single-  and  multiple-pulse  exposures,  for  absorption  coefficients  corresponding  to 
10.6  pm  and  2.01  pm  radiation,  and  for  exposure  durations  ranging  from  1.7  ns  to 
5  sec.  In  all  cases  the  results  agreed  with  the  previous  calculations  to  the  number 
of  decimal  figures  printed  in  the  original  mainframe  output  (usually  6). 

After  re-installing  the  experimental  set-up,  we  first  made  a  few  single¬ 
pulse  exposures  to  reconfirm  our  earlier  results  and  to  refamilarize  ourselves 
with  the  experimental  protocol.  Using  a  TEMqo  beam  with  a  1/e  diameter  of  3.32 
mm  we  exposed  corneas  with  peak  energy  densities  of  442,  397  and  329  mJ/cm^. 
The  first  two  exposures  caused  superficial  lesions  and  the  third  exposure  did  not. 
In  our  earlier  study^- 1  we  had  obtained  a  threshold  of  360  mJ/cm^;  thus,  these 
exposures  confirmed  that  the  apparatus  and  our  methodology  give  results  for 
single  pulses  that  are  consistent  with  those  obtained  previously. 

We  set  out  to  refine  the  threshold  for  2  pulses  at  10  Hz  but  encountered 
difficulty  obtaining  consistent  pulses.  The  second  pulse  either  had  an  amplitude 
that  was  significantly  lower  than  the  first,  or  was  missing  altogether.  Initial 
attempts  to  solve  this  problem  were  unsuccessful.  A  major  difficulty  was  that 
Boston  Laser  has  been  out  of  business  for  several  years  and  the  documentation  for 
the  laser  is  very  sparse  ( it  lacks  most  circuit  and  gas  flow  diagrams).  We  did  find 
that  we  were  able  to  obtain  consistent  pulses  up  to  5  Hz.  Thus  we  temporarily 
proceeded  with  determining  the  damage  threshold  for  a  sequence  of  2  pulses  at  5 
Hz.  The  result,  275  mJ/  cmVpulse,  lies  between  our  earlier  results  for  2  pulses  at  1 
Hz  (300  mJ/  cmVpulse)  and  2  pulses  at  10  Hz  (200  mJ/  cmVpulse)^’  7.  This 
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threshold  exposure  resulted  in  a  calculated  peak  temperature  increase  on  the 
beam  axis  10  pm  beneath  the  anterior  tear  layer  of  26.2  C  (see  Fig.  2). 
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Figure  2:  Temperature  history  for  threshold  exposure  (2  pulses  at  5  Hz). 

After  completing  this  determination  we  returned  to  solving  the  problem  of 
inconsistent  pulsing  at  frequencies  above  5  Hz.  We  cleaned  and  refurbished  the 
spark  gap  and  verified  that  the  capacitors  were  charging  and  discharging 
properly.  Most  of  this  work  was  done  subsequent  to  the  closing  date  for  this  report; 
however  the  problem,  which  we  found  involved  gas  flow,  is  solved.  In  order  to 
verify  operation  at  10  Hz  we  used  the  Scientech  detector  to  measure  the  energy  of 
ten  single  pulses  and  the  integrated  energy  of  ten  2-pulse  sequences.  We  obtained 
an  average  value  (±S.D.)  of  29.9Q±1.79  mJ  for  the  single  pulses  and  30.75±1.06  m  J 
for  one-half  the  integrated  energy  of  the  2-pulse  sequences.  These  values  are 
statistically  indistingmshable  (t  =  0.11).  We  are  now  in  the  process  of  completing 
the  damage  determinations  given  in  the  Statement  of  Work. 

Our  Contracting  Officer’s  Representative,  Mr.  Bruce  Stuck,  and  Dr.  David 
Sliney  of  the  United  States  Army  Center  for  Health  Promotion  and  Preventive 
Medicine  visited  APL  on  July  31  to  discuss  progress  and  future  plans.  We 
discussed  several  factors  to  consider  in  planning  future  research.  Among  these 
were  the  importance  of  gathering  data  in  the  1.3  to  1.5  pm  spectral  region  where 
the  radiation  is  more  penetrating.  In  this  spectral  region  particularly,  and 
perhaps  in  the  region  near  2  pm,  it  was  emphasized  that  it  would  be  important  to 
consider  late  effects  as  well  as  immediate  effects  of  exposures.  As  we  progress 
with  the  threshold  determinations  we  will  confer  again  with  Mr.  Stuck  to 
formulate  the  next  steps. 

Conclusions 

The  program  we  developed  for  calculating  temperature  histories  resulting 
from  exposure  to  a  laser  beam  having  a  Gaussian  irradiance  profile  was  updated 
to  run  on  a  desktop  computer.  The  revision  required  new  routines  for  computing 
complementary  error  functions  and  for  evaluating  the  definite  integral.  The 
program  was  tested  extensively  and  it  was  shown  to  give  results  in  exact 
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agreement  with  the  earlier  version.  We  reconfirmed  our  earlier  damage 
threshold  for  a  single  80  ns  pulse  and  determined  the  damage  threshold  for  a 
sequence  of  two  80  ns  pulses  at  a  pulse  repetition  frequency  of  5  Hz.  The  result 
was  consistent  with  our  earlier  determinations  of  2  pulses  at  1  Hz  and  2  pulses  at 
10  Hz.  Problems  of  achieving  uniform  pulse  amplitudes  at  10  Hz  were  solved. 
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Appendix  2 


PROGRAM  CALCULATES  TEMPERATURE-TIME  fflSTORIES  FOR  RADIATION 
INCIDENT  ON  A  SEMI-INFINITE  ABSORBING  SLAB 
RADIATION  CAN  BE  EITHER  CW  OR  SEQUENCES  OF  PULSES 
RADIATION  AND  CONVECTIVE  HEAT  TRANSFER  ARE  IGNORED 
PROGRAM  ADAPTED  BY  R.  MCCALLY  (FEBRUARY  14-MARCH  18  1996)  FROM 
MAINFRAME  VERSION  WRITTEN  BY  S.  FAVIN 

SO  =  IRRADIANCE  (W/(CM*CM)) 

KAPPA  =  THERMAL  CONDUCTIVITY  (CAL/(CM*CM*DEGC)) 

R  =  RADIAL  DIST  FROM  BEAM  AXIS  (CM) 

Z  =  DEPTH  (CM) 

RE  =  1/E  BEAM  RADIUS  (CM) 

NPTS  =  NUMBER  OF  CALC  POINTS 
TAU  =  PULSE  DURATION  (SEC) 

C  =  HEAT  CAPACITY  (CAL/(GM*DEGC)) 

GAM  =  INVERSE  ABSORBTION  LENGTH  ( 1/CM) 

DELT  =  1/PRF  (SEC)  —  DELT  MUST  BE  >=  TAU 
NP  =  NUMBER  OF  PULSES 

ND  =  NUMBER  OF  INTERVALS  CALCULATED  DURING  EACH  PULSE 
NRUNS  =  NUMBER  OF  SETS  OF  INPUT  DATA 

IRULE  -  Choice  of  quadrature  rule  in  DQDAG.  (Input) 

A  Gauss-Kronrod  rule  is  used  with 
7-15  points  if  IRULE  =  1 
10-21  points  if  IRULE  =  2 
15-31  points  if  IRULE  =  3 
20-41  points  if  IRULE  =  4 
25-51  points  if  IRULE  =  5 
30-61  points  if  IRULE  =  6 
IRULE  =  2  is  recommended  for  most  functions. 

If  the  function  has  a  peak  singularity  use  IRULE  =  1. 

If  the  function  is  oscillatory  use  IRULE  =  6. 

PROGRAM  IRTEMP 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  KAPPA 

COMMON  Z0,0NE,TW0,TEN,HUN,SPI,W,ZA,GTA,AZ,E2,R2,  XO,  TEGZ,  G2FA 
DIMENSION  TIME(8001),TEMP(8001),TEMPN(8001) 

DIMENSION  IP(IOOO) 

CHARACTER*!  TAB 
EXTERNAL  FF 

INN  =  7  !  DEFINES  INPUT  LOGICAL  UNIT 

lOUT  =  8  !  DEFINES  OUTPUT  LOGICAL  UNIT  FOR  FINAL  DATA  FILE 

ISCREEN  =  6  !  DEFINES  LOGICAL  UNIT  FOR  OUTPUT  TO  SCREEN 

ICOUNT  =  0  IINITIALIZE  COUNTER  FOR  TIMES  INPUT  FILE  IS  ACCESSED 

DEFINE  CONSTANTS 
ZO  =  O.ODO 
ONE=  l.ODO 
TWO  =  2.0D0 
TEN  =  lO.ODO 
HUN  =  lOO.ODO 


12 


on  n  nn^n  nn 


11/26/96  5:12  PM 


PI  =  4.0D0  *  DATAN(  ONE ) 

SPI  =  DSQRT(  PI ) 

TAB  =  CHAR(9)  !  TAB  IS  ASC  9 

OPEN  I/O  FILES 

OPEN  (UN1T=INN,  FILE  ='IRINPUT’,  STATUS='OLD’,ACCESS='SEQUENTIAL',  BLANK='NULL’) 
OPEN  (UNIT=IOUT,  FILE  ='IROUT’,  STATUS='UNKNOWN',ACCESS='SEQUENTIAL') 

m  CONTINUE 

GET  INPUT  DATA 

CALL  INPUT(S0, KAPPA, R,Z,RE,NPTS,TAU,C,GAM,DELT,NP,ND,RHO,NRUNS,ICOUNT, INN, 
1I0UT,ISCREEN,AERR,RERR,IRULE) 

A2  =  RHO*C/(4.0D0*KAPPA) 

TSIG  =  RE*RE 
R2  =  R*R 

FRONT  =  SO*GAM/(8.36DO*RHO*C) 

ZA  =  Z**2*  A2 
AZ  =  DSQRT(ZA) 

TA  =  TWO  *  DSQRT(  A2 ) 

GT  A  =  GAM/T  A 

XO  =  TWO*A2*Z/GAM 

TEGZ  =  TWO  *  DEXP(-GAM*Z) 

G2FA  =  GAM**2  /  (4.D0*A2) 

W  =  A2  *  TSIG 

CALCULATE  TIME  POINTS 
TIME(l)  =  ZO 
DO30I=l,ND 
TIME(I+1)  =  FTAU/ND 
30  CONTINUE 
NDl  =  ND  +  1 
N2  =  2*ND 

C  FOR  SINGLE  PULSES  THIS  LOOP  GENERATES  ND  TIME  POINTS  BETWEEN  TAU  AND  DELT 
C  IF  NPTS  <2*ND  IT  GENERATES  NPTS-ND  POINTS  WITH  SPACING  (DELT-TAU)/ND 
DO  32  I=ND1,N2 

TIME(I+1)  =  TAU  +  (I-ND)*(DELT  -  TAU)/ND 
32  CONTINUE 
ND2  =  N2  +  1 

C  FOR  SINGLE  PULSES  THIS  LOOP  GENERATES  NPTS-2*ND  TIME  POINTS  AFTER  DELT  (SPACI 
C  SAME  AS  DURING  PULSE) 

DO  36  I=ND2,NPTS  !  DOES  NOT  EXECUTE  IF  NPTS<2*ND+ 1 

D0  34J=1,NP 
K  =  J*N2 

IF(I  .LE.  K)  GO  TO  34 
TIME(I+1)  =  J*DELT  +  TIME(I+1-K) 

34  CONTINUE 
36  CONTINUE 
TEMP(l)  =  ZO 
WRrrE(ISCREEN,4) 

WRITE(IOUT,5)TAB,TAB,TAB 

C 

C  BEGIN  TEMPERATURE  CALCULATION 
DO  100 1=2, NPTS 
TL  =  TIME(1) 

IF(TIME(I)-TAU  .GT.  TL)  TL=TIME(I)-TAU 
TU  =  TIME(I) 
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CALL  DQDAG  (FF,  TL,  TU,  AERR,  RERR,  IRULE,  PSI,  ERREST) 

TEMP(I)  =  FRONT  *  PSI 

WRITE(ISCREEN,6)  I,PSI,TIME(I),TEMP(I), ERREST 
WRITE(IOUT,7)  I,TAB,PSLTAB,TIME(I),TAB,TEMP(I), TAB, ERREST 
100  CONTINUE 
C 

C  HANDLE  MULTIPLE  PULSES 
IF(NP.EQ.l)  GO  TO  400 
IP(1)=1 
DO  50  J=2,NP 
IP(J)  =  (J-1)  *  ND  *  2 
50  CONTINUE 
DO  300  I=1,NPTS 
TEMPN(I)  =  TEMP(I) 

IF(TIME(I)  .LT.  DELT)  GO  TO  300 
D0  200J=2,NP 
IF(I  .LE.  IP(J))  GO  TO  200 
L  =  I  -  IP(J) 

TEMPN(I)  =  TEMPN(I)  +  TEMP(L) 

200  CONTINUE 
300  CONTINUE 
GO  TO  375 
350  DO  375  I=1,NPTS 
TEMPN(I)  =  TEMP(I) 

375  CONTINUE 
C 

C  WRITE  TO  FINAL  DATA  FILE  (ONLY  IF  MULTIPLE  PULSES) 

WRrrE(IOUT,380)TAB,TAB  !  WRITE  HEADER  FOR  RESULTS 

DO  400  I=1,NPTS 

WRn'E(IOUT,390)  I,TAB,TIME(I),TAB,TEMPN(I) 

400  CONTINUE 

IF  (ICOUNT  .LT.  NRUNS)  GO  TO  1000 
99  STOP 

C  3  FORMATC  A2  =  ',1PE14.5,'  TSIG  =  ’,1PE14.5,’  R2  =  ’,1PE14.5,'  FRONT  =  MPE14.5) 

4  FORMAT(T3,T,Tl  1, 'PSI, T25,TIME’,T36,'TEMPERATURE',T52, 'ERROR') 

5  FORMAT('r,Al,'PSr,Al,'TIME',Al, 'TEMPERATURE', Al, 'ERROR') 
6FORMAT(I5,lPG14.6,lPG14.6,lPG14.6,lPG14.6) 
7FORMAT(I5,Al,lPG14.7,Al,lPG14.7,Al,lPG14.7,Al,lPG14.7) 

380FORMAT('r,Al,'TIME',Al,'TEMPERATURE') 

390FORMAT(I5,A1,1PE14.8,A1,1PE14.8) 

END 

C 

C 

FUNCTION  FF(X) 

C  THIS  DESCRIBES  THE  INTEGRAND 
IMPLICIT  REAL*8(A-Z) 

COMMON  ZO,ONE,TWO,TEN,HUN,SPI,W,ZA,GTA,AZ,E2,R2,  XO,  TEGZ,  G2FA 
C 

IF(X  .EQ.  ZO)  GO  TO  50 
C 

X2  =  DSQRT(  X ) 

XW  =  ONE  +  XAV 
EE  =  -R2/XW 
El  =  ZO 

IF(DABS(EE)  .LE.  HUN)  El  =  DEXP(  EE ) 

C 

E2  =  -ZA/X 
C 


14 


11/26/96  5:12  PM 


Cl  =GTA*X2 
C2  =  AZ/X2 
ARGl  =  Cl  +  C2 
ARG2  =  C1  -C2 
HI  =  HFUN(  ARGl ) 

H2  =  HFUN(  ARG2 ) 

C 

IF(X  .GT.  XO)  GO  TO  10 

C...  FOR  X  LESS  THEN  XO  AND  NOT  ZERO  6/2/80  THE  REVISION 
FO  =  El/XW 

G1  =  TEGZ  *  DEXP(G2FA*X) 

G2=ZO 

IF(DABS(E2)  .GT.  HUN)  GO  TO  30 
HI  =  HFUN(  DABS(ARGl) ) 

H2  =  HFUN(  DABS(ARG2) ) 

G2  =  (H1-H2) 

C 

30  FF  =  FO*(Gl  +  G2) 

RETURN 

C 

10  CONTINUE 
FF  =  E1  *(H1+H2)/XW 
RETURN 

C...  FOR  X=  0  ONLY 
50  CONTINUE 
FF  =  TEGZ  *  DEXP(-R2) 

RETURN 

END 

C 

C 

FUNCTION  HFUN  ( ARG ) 

IMPLICIT  REAL*8(A-Z) 

COMMON  Z0,0NE,TW0,TEN,HUN,SPI,W,ZA,GTA,AZ,E2,R2,  XO,  TEGZ,  G2FA 
DATAHALF,ONEH/0.5D0,  L50D0/ 

C 

Y  =  DABS(  ARG ) 

Y2  =  Y*Y 
HFUN  =  ZO 
BB  =  Y2  +  E2 

IF(ARG  .GE.  -13.0D0)  GO  TO  10 
Bl=ZO 

IF(DABS(BB)  .LE.  174.0D0)  B1  =  TWO  *  DEXP(  BB  ) 

B2  =  ZO 

IF(DABS(E2)  .LT.  174.0D0)  B2  =  DEXP(E2)  / 

*  (SPI*(Y+HALF/(Y+ONE/(Y+ONEH/(Y+TWO/Y)))) ) 

HFUN  =  B1-B2 
RETURN 
10  CONTINUE 
IF(ARG  .GE.  ZO)  GO  TO  20 

IF(DABS(BB)  .LT.  174.0D0)  HFUN  =  DEXP(BB)*(TWO-DERFC(Y)) 

RETURN 
20  CONTINUE 
IF(ARG  .GT.  13.0)  GO  TO  30 

IF(DABS(BB)  .LT.  174.0D0)  HFUN  =  DEXP(BB)*DERFC(ARG) 

RETURN 
30  CONTINUE 

IF(DABS(E2)  .GT.  174.0D0)  GO  TO  40 
Z  =  ARG 
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HFUN  =  SPI*(Z+HALF/(Z+ONE/(Z+ONEH/(Z+TWO/Z)))) 
IF(HFUN  .NE.  ZO)  HFUN  =  DEXP(E2)/HFUN 
40  RETURN 
END 
C 
C 

FUNCTION  DERFC(X) 

C 

C  RETURNS  THE  COMPLIMENTARY  ERROR  FUNCTION 


c 

c 

c 

c 

c 


checks  with  Mathematica  to  15  significant  figures 
in  range  (.001  <=  x  <=  4.0)  (Maximum  printout  from  Mathematica 
and  maximum  range  checked  3/18/96  ^M) 

implicit  real*  8  (a-h,o-z) 

DATA  Z0,HALF,ONE/0.0D0,0.5D0, 1  .ODO/ 

IF(X.LT.zO)THEN 

DERFC=ONE+GAMMP(half,X**2) 

ELSE 

DERFC=GAMMQ(half,X**2) 

ENDIF 

RETURN 

END 


FUNCTION  GAMMP(A,X) 
implicit  real*  8  (a-h,o-z) 

DATA  Z0,HALF,ONE/0.0D0,0.5D0,  l.ODO/ 
IF(X.LT.zO.OR.A.LE.zO)PAUSE 
IF(X.LT.A+one)THEN 
CALL  GSER(GAMSER,A,X,GLN) 
GAMMP=GAMSER 
ELSE 

CALL  GCF(GAMMCF,A,X,GLN) 
GAMMP=one-GAMMCF 
ENDIF 
RETURN 
END 


FUNCTION  GAMMQ(A,X) 
implicit  real*  8  (a-h,o-z) 

DATA  Z0,HALF,ONE/0.0D0,0.5D0,  l.ODO/ 
IF(X.LT.ZO.OR.A,LE.ZO)PAUSE 
IF(X.LT.A+ONE)THEN 
CALL  GSER(GAMSER,A,X,GLN) 
GAMMQ=ONE-GAMSER 
ELSE 

CALL  GCF(GAMMCF,A,X,GLN) 
GAMMQ=GAMMCF 
ENDIF 
RETURN 
END 


SUBROUTINE  GSER(GAMSER,A,X,GLN) 
implicit  real*8  (a-h,o-z) 

DATA  Z0,HALF,ONE/0.0D0,0.5D0,  l.ODO/ 

PARAMETER  (ITMAX=400,EPS=10D-17) 

gin  =  0.5723649429247001d0  !  ln(sqrt(pi))  =  gamma(l/2) 
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IF(X.LE.zO)THEN 
IF(X.LT.zO)PAUSE 
GAMSER=zO 
RETURN 
END  IF 
AP=A 

SUM=one/A 
DEL=SUM 
DO  11  N=1,ITMAX 
AP=AP+one 
DEL=DEL*X/AP 
SUM=SUM+DEL 

IF(dABS(DEL).LT.dABS(SUM)*EPS)GO  TO  1 
11  CONTINUE 

PAUSE  'A  too  large,  FTMAX  too  small' 

I  GAMSER=SUM*dEXP(-X+A*dLOG(X)-GLN) 

RETURN 

END 

c 

SUBROUTINE  GCF(GAMMCF,A,X,GLN) 
implicit  real*  8  (a-h,o-z) 

DATA  Z0,HALF,ONE/0.0D0,0.5D0,1.0D0/ 

PARAMETER  (ITMAX=400,EPS=  lOd- 17) 

gin  =  0.5723649429247001d0  !  ln(sqrt(pi))  =  gamma(l/2) 

GOLD=zO 

aO=one 

A1=X 

BO=zO 

Bl=one 

FAC=one 

DO  11  N=l,rTMAX 
AN=dFLOAT(N) 

ANA=AN-A 

A0=(A1+A0*ANA)*FAC 
B0=(B  1+B0*ANA)*FAC 
ANF=AN*FAC 
A1=X*A0+ANF*A1 
B1=X*B0+ANF*B1 
IF(Al.NE.zO)THEN 
FAC=one/Al 
G=B1*FAC 

IF(DABS((G-GOLD)/G).LT.EPS)GOTO  1 
GOLD=G 
ENDIF 

II  CONTINUE 

PAUSE  'A  too  large,  FTMAX  too  small' 

1  GAMMCF=dEXP(-X+A*dLOG(X)-GLN)*G 
RETURN 
END 


SUBROUTINE  INPUT(SO,KAPPA,R,Z,RE,NPTS,TAU,C,GAM,DELT,NP,ND,RHO,NRUNS, 
lICOUNT,INN,IOUT,ISCREEN,AERR,RERR,IRULE) 

C - 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  KAPPA 

C  INTEGER*4  NPTS,NP,ND,NRUNS,ICOUNT,INN,IOUT,ISCREEN,IRULE,I 
CHARACrER*9  DAY 
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CHARACTER*  11  ADI 
CHARACTER*30  AIN2 


PURPOSE;  READ  INPUT  FILE  IRINPUT  AND  WRITE  PARAMETERS. 


INPUT;  INN,IOUT,ISCREEN 


OUTPUT;  SO,KAPPA,R,Z,RE,NPTS,TAU,C,GAM,DELT,NP,ND, 
RHO,NRUNS,ICOUNT,AERR,RERR,IRULE. 

ICOUNT=ICOUNT+l  !  INCREMENT  FILE  ACCESS  COUNTER 

READ  FILE  INN. 

READ(INN;(A30)’,END=99,IOSTAT=I)AIN2  !  READ  HEADER 

READ(INN,'(A30)',END=99,IOSTAT=I)AIN2  !  READ  — 

READ(INN,'(A1 1)',END=99,I0STAT=I)AD1  !  READ123  ETC 

READ(INN;(  A 1 1 ,15)',END=99,IOSTAT=I)  AD  1  ,NRUNS 
READ(INN;(A  1 1 , 1PG14.6)',END=99,I0STAT=I)  ADI  ,S0 
READ(INN,'(  A1 1 , 1PG14.6)',END=99,I0STAT=I)  AD  1  ,TAU 
READ(INN,'(  A 1 1  ,F  14.9)',END=99,IOSTAT=I)  AD  1  ,RE 
READ(INN,'(A11,F14.9)’,END=99,I0STAT=I)AD1,Z 
READ(INN;(  A 1 1  ,F  1 4.9)’,END=99,IOSTAT=I)  AD  1  ,R 
READ(INN,'(A  1 1 ,15)',END=99,IOSTAT=I)  AD  1  ,NP 
READ(INN,'(A1 1 , 1PG14.6)',END=99,I0STAT=I)AD  1  ,DELT 
READ(INN,'(  A 1 1 ,15)’,END=99,IOSTAT=I)  AD  1  ,ND 
READ(INN,’(A11,I5)’,END=99,I0STAT=I)AD1,NPTS 
READ(INN,'(  A 1 1  ,F  14.9)’,END=99,IOSTAT=I)  AD  1  ,GAM 
READ(INN,'(  A 1 1  ,F  14.9)’,END=99,IOSTAT=I)  AD  1  ,C 
READ(INN,'(A1 1,F14.9)',END=99,I0STAT=I)AD1, KAPPA 
READ(INN,'(  A 1 1  ,F  14.9)',END=99,IOSTAT=I)  AD  1  ,RHO 

. INTEGRATION  PARAMETERS; 

READ(INN,'(A1 1,F14.9)‘,END=99,I0STAT=I)AD1,AERR 
READ(INN;(  A1 1,  IPD  14.9)’,END=99,IOSTAT=I)AD  1  ,RERR 
READ(INN,’(A11,I5)',END=99,I0STAT=I)AD1,IRULE 

WRITE  INPUT  TO  FILE  IROUT  ( DISK  DATA  FILE) 

WRITE(IOUT,100) 

WRITE(IOUT,105) 

CALLDATE(DAY) 

WRITE(IOUT,  1 10)DAY 
AIN2  =  ’  INPUT  PARAMETERS’ 

WRrTE(IOUT,110)AIN2 

WRITE(IOUT,150)NRUNS 

WRITE(IOUT,151)SO 

WRITE(IOUT,  152)T  AU 

WRITE(IOUT,153)RE 

WRITE(IOUT,154)Z 

WRITE(IOUT,166)R 

WRITE(IOUT,155)NP 

WRrrE(IOUT,156)DELT 

WRITE(IOUT,157)ND 

WRrrE(IOUT,158)NPTS 

WRITE(IOUT,159)GAM 

WRITE(IOUT,160)C 

WRrTE(IOUT,161)KAPPA 

WRITE(IOUT,  162)RHO 
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AIN2  =  ' INTEGRATION  PARAMETERS:  ’ 
WRITE(IOUT,110)AIN2 
WRITE(IOUT,  163)IRULE 
WRITE(IOUT,164)AERR 
WRITE(IOUT,  165)RERR 
C 

C  WRITE  INPUT  TO  SCREEN 
C 

WRITE(ISCREEN,100) 

WRITE(ISCREEN,  1 10)DAY 
AIN2  =  ’  INPUT  PARAMETERS' 

WRITE(ISCREEN,  1 10)  AIN2 
WRITE(ISCREEN,  150)NRUNS 
WRITE(ISCREEN,15 1)S0 
WRITE(ISCREEN,  152)TAU 
WRITE(ISCREEN,153)RE 
WRITE(ISCREEN,  154)Z 
WRITE(ISCREEN,  166)R 
WRITE(ISCREEN,155)NP 
WRITE(ISCREEN,  156)DELT 
WRITE(ISCREEN,  157)ND 
WRITE(ISCREEN,  158)NPTS 
WRITE(ISCREEN,  159)GAM 
WRITE(ISCREEN,160)C 
WRITE(ISCREEN,161)KAPPA 
WRITE(ISCREEN,  162)RHO 
AIN2  =  ' INTEGRATION  PARAMETER:  ' 
WRITE(ISCREEN,1 10)AIN2 
WRITE(ISCREEN,  163)IRULE 
WRITE(ISCREEN,  164)AERR 
WRITE(ISCREEN,  165)RERR 

99  CONTINUE 
C 

C  FINAL  WRITE  STATEMENTS 
C 

WRrrE(IOUT,120) 

WRITE(ISCREEN,120) 

RETURN 

C 

100  FORMAT(/'  —WE  ARE  IN  IRINPUT  —  '/) 

105  FORMAT(/'  —  FILE  IROUT—  ’/) 

no  FORMAT(A35) 

120  FORMAT(/'  —LEAVING  IRINPUT—  ’/) 

130  FORMATC ... ’) 

150  FORMATC  NRUNS  =',I5) 

151  FORMATC  SO  = ’,1PG14.5,' WATTS/CM**2’) 

152  FORMATC  TAU  =  ',1PG14.5,’  SEC) 

153  FORMATC  RE  =  ’,F14.5;  CM') 

154  FORMATC  Z  =  ',F14.5;  CM') 

155  FORMATC  NP  =',I5,'  PULSES') 

156  FORMATC  DELT  = ',1PG14.5,’ SEC) 

157  FORMATC  ND  = ',15,'  INTERVALS/PULSE’) 

158  FORMATC  NPTS  =  ',15) 

159  FORMATC  GAM  =  ’,F14.5,'  1/CM') 

160  FORMATC  C  =  ',F14.5,'  CAL/(GM*DEGC)') 

161  FORMATC  KAPPA  = ',F14.5,' CAL/(CM**2*DEGC)') 

162  FORMATC  RHO  =  ',F14.5,’  GM/CM**2’) 

163  FORMATC  IRULE  =',I5) 
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164  FORMATCAERR  =',F14.5) 

165  FORMATCRERR  =’,1PD14.5) 

166  FORMATCR  = ’,F14.5,' CM') 

END 

C 

c - - - 

C  IMSL  Name:  QDAG/DQDAG  (Single/Double  precision  version) 

C 

C  Computer:  PCDSMS/DOUBLE 
C 

C  Revised:  January  29,  1985 
C 

C  Purpose:  Integrate  a  function  using  a  globally  adaptive 
C  scheme  based  on  Gauss-Kronrod  rules. 

C 

C  Usage:  CALL  QDAG  (F,  A,  B,  ERRABS,  ERRREL,  IRULE,  RESULT, 
C  ERREST) 

C 

C  Arguments: 

C  F  -  User-supplied  FUNCTION  to  be  integrated.  The  form  is 
C  F(X),  where 

C  X  -  Independent  variable.  (Input) 

C  F  -  The  function  value.  (Output) 

C  F  must  be  declared  EXTERNAL  in  the  calling  program. 

C  A  -  Lower  limit  of  integration.  (Input) 

C  B  -  Upper  limit  of  integration.  (Input) 

C  ERRABS  -  Absolute  accuracy  desired.  (Input) 

C  ERRREL  -  Relative  accuracy  desired.  (Input) 

C  IRULE  -  Choice  of  quadrature  rule.  (Input) 

C  A  Gauss-Kronrod  rule  is  used  with 

C  7-15  points  if  IRULE  =  1 

C  10-21  points  if  IRULE  =  2 

C  15-31  points  if  IRULE  =  3 

C  20-41  points  if  IRULE  =  4 

C  25-51  points  if  IRULE  =  5 

C  30-61  points  if  IRULE  =  6 

C  IRULE  =  2  is  recommended  for  most  functions. 

C  If  the  function  has  a  peak  singularity  use  IRULE  =  1 . 

C  If  the  function  is  oscillatory  use  IRULE  =  6. 

C  RESULT  -  Estimate  of  the  integral  from  A  to  B  of  F.  (Output) 

C  ERREST  -  Estimate  of  the  absolute  value  of  the  error.  (Output) 

C 

C  Remarks: 

C  1.  Automatic  workspace  usage  is 
C  QDAG  2500  units,  or 

C  DQDAG  4500  units. 

C  Workspace  may  be  explicitly  provided,  if  desired,  by  use  of 
C  Q2AG^Q2AG.  The  reference  is 

C  CALL  Q2AG  (F,  A,  B,  ERRABS,  ERRREL,  IRULE,  RESULT, 

C  ERREST,  MAXSUB,  NEVAL,  NSUBIN,  ALIST,  BLIST, 

C  RLIST,  ELIST,  lORD) 

C  The  additional  arguments  are  as  follows: 

C  MAXSUB  -  Number  of  subintervals  allowed.  (Input) 

C  A  value  of  500  is  used  by  QDAG. 

C  NEVAL  -  Number  of  evaluations  of  F.  (Output) 

C  NSUBIN  -  Number  of  subintervals  generated.  (Output) 

C  ALIST  -  Array  of  length  MAXSUB  containing  a  list  of  the  NSUBIN 
C  left  endpoints.  (Output) 
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BUST  -  Array  of  length  MAXSUB  containing  a  list  of  the  NSUBIN 
right  endpoints.  (Output) 

RLIST  -  Array  of  length  MAXSUB  containing  approximations  to 
the  NSUBIN  integrals  over  the  intervals  defined  by  ALIST, 

BUST.  (Output) 

BUST  -  Array  of  length  MAXSUB  containing  the  error  estimates  of 
the  NSUBIN  values  in  RLIST.  (Output) 
lORD  -Array  of  length  MAXSUB.  (Output) 

Let  K  be 

NSUBIN  if  NSUBIN  .LE.  (MAXSUB/2+2) 

MAXSUB+ 1  -NSUBIN  otherwise. 

The  first  K  locations  contain  pointers  to  the  error 
estimates  over  the  corresponding  subintervals,  such  that 
ELIST(IORD(l)), ...,  ELIST(IORD(K))  form  a  decreasing 
sequence. 

2.  Informational  errors 
Type  Code 

4  1  The  maximum  number  of  subintervals  allowed  has  been 
reached. 

3  2  Roundoff  error,  preventing  the  requested  tolerance  from 
being  achieved,  has  been  detected. 

3  3  A  degradation  in  precision  has  been  detected. 

3.  If  EXACT  is  the  exact  value,  QDAG  attempts  to  find  RESULT 

such  that  ABS(EXACT-RESULT)  .LE.  MAX(ERRABS,ERRREL*ABS(EXACT)). 
To  specify  only  a  relative  error,  set  ERRABS  to  zero.  Similarly, 
to  specify  only  an  absolute  error,  set  ERRREL  to  zero. 

GAMS:  H2alal 

Chapter:  MATH/LIBRARY  Integration  and  Differentiation 

Copyright:  1985  by  IMSL,  Inc.  All  Rights  Reserved. 

Warranty:  IMSL  warrants  only  that  IMSL  testing  has  been  applied 
to  this  code.  No  other  warranty,  expressed  or  implied, 
is  applicable. 
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